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ABSTRACT 


We  present  progress  in  the  development  of  a  new  approach  to  develop  and  evaluate  earth  models  at  the  regional 
scale  that  utilizes  full  waveform  seismograms. 

We  have  recently  implemented  a  non-linear  three-dimensional  (3-D)  Born  approximation  and  constructed  a  3-D 
shear  velocity  model  of  a  subregion  of  Southeast  Asia  (longitude  75  to  150  degrees  and  latitude  0  to  45  degrees) 
based  on  this  approach  for  the  forward  modeling  part  of  the  problem  and  linear  Born  kernels  for  the  inverse  part. 

Our  initial  model  in  the  large  region  (longitude  30  to  150  and  latitude  -10  to  60  degrees)  is  derived  from  a  large  data 
set  of  teleseismic  surface  waveforms  (fundamental  mode  and  overtones  in  the  period  range  60-300  s)  using  the 
Nonlinear  Asymptotic  Coupling  Theory  (NACT),  an  approach  used  successfully  for  global  and  regional  mantle 
tomography  at  Berkeley  since  1995.  In  the  subregion  of  study,  our  “N-Bom”  model  is  parameterized  at  relatively 
short  wavelengths  on  the  order  of  200  km. 

The  NACT  approach  assumes  two-dimensional  (2-D)  sensitivity  kernels,  confined  to  the  vertical  plane  containing 
the  source  and  the  receiver,  and  is  adequate  for  the  development  of  a  smooth  velocity  model.  We  first  developed  a 
model  for  Southeast  Asia  using  NACT.  Using  this  model  as  a  starting  model,  we  performed  one  iteration  using  the 
N-Bom  approach.  This  approach  computes  synthetic  seismograms  including  the  effects  of  single  scattering  (linear) 
as  well  as  a  non-linear  part,  as  in  the  standard  “path  average  approximation,”  that  accounts  for  large  accumulated 
phase  delays  on  paths  that  sample  large  scale  smooth  anomalies.  The  inversion  kernels  are  complete  3-D  Bom 
kernels. 

A  regional  version  of  the  Spectral  Element  Method  (SEM)  code,  RegSEM.  1 ,  has  been  completed.  This  version  of 
the  code  accepts  a  non-conformal  grid,  uses  PML  (Perfectly  Matched  Layers)  at  the  borders  of  the  region,  and 
includes  general  3-D  anisotropy,  Moho  and  surface  topography,  ocean  bathymetry,  attenuation,  and  ellipticity.  We 
are  computing  synthetics  in  our  N-Bom  model  and  will  perform  one  iteration  of  inversion  using  3-D  Born  sensitivity 
kernels  to  produce  the  next  generation  3-D  model  of  upper  mantle  stmcture  in  the  subregion  defined  above. 

In  order  to  obtain  finer  velocity  stmctures  beneath  the  region  of  our  study,  we  are  performing  forward  modeling  at 
shorter  periods  using  the  method  of  frequency-wave  number  integration  (FKI).  Using  our  N-Bom  3-D  model,  we 
add  sedimentary  layers  in  the  cmst  and  change  their  thickness  to  improve  the  fit  to  the  observed  seismograms,  and 
then,  based  on  the  refined  velocity  stmctures,  we  conduct  our  own  moment  tensor  inversions,  as  intermediate  steps. 
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OBJECTIVES 


The  primary  objective  of  this  research  is  to  develop  and  apply  an  approach  to  utilize  increasingly  advanced 
theoretical  frameworks  and  numerical  methods  to  obtain  improved  regional  seismic  structure  calibration. 
Specifically,  a  large-scale  regional  Eurasian  model  has  been  developed  from  a  large  data  set  of  seismic  waveforms 
using  the  path-average  approximation  (PAVA)  and  NACT  (Li  and  Romanowicz,  1995),  which  are  well-developed 
normal-mode  based  approaches  that  consider  one-dimensional  (1-D)  and  2-D  waveform  sensitivity  respectively 
along  the  great-circle  path  between  source  and  receiver.  We  refine  this  model  in  a  smaller  region  using  a  linear 
implementation  of  Bom  single-scattering  theory  (Capdeville,  2005),  which  more  accurately  represents  the  3-D 
sensitivity  of  the  seismic  wavefield,  for  the  inverse  part  of  the  problem  and  using  a  non-linear  modified  3-D  Bom 
approximation  for  the  forward  part  of  the  problem  (Panning  et  al.,  2006).  Additionally,  we  have  developed  and 
implemented  a  regional  version  of  SEM,  a  numerical  approach  that  accurately  models  wave  propagation  in  a 
complex  3-D  Earth  (e.g.,  Faccioli  et  al.,  1996;  Komatitsch  and  Vilotte,  1998;  Komatitsch  and  Tromp,  1999).  Finally, 
we  plan  to  utilize  a  novel  approach  with  stacked  sources  (Capdeville  et  al.,  2003)  to  further  speed  computation. 

A  further  objective  of  this  research  is  to  perform  validation  and  improved  calibration  of  the  model  described  above 
using  a  variety  of  approaches  and  data  sets,  including  ground-tmth  data  sets  from  the  Knowledge  Base.  Specifically, 
we  plan  to  apply  teleseismic  receiver  function  modeling,  regional  broadband  data  forward  modeling,  and  surface 
wave  group  velocity  measurements  to  test  and  improve  the  model  using  data  not  included  in  the  original  inversion. 
While  these  additional  data  sets  can  help  improve  all  aspects  of  the  model,  we  anticipate  the  greatest  improvement  in 
the  shallow  stmcture,  particularly  the  cmst,  which  is  not  as  well  constrained  by  the  longer-period  data  used  in  the 
initial  model  development. 

This  research  can  then  serve  as  a  proof  of  concept  for  applying  a  similar  approach  to  the  calibration  of  seismic 
structure  in  other  regions  of  Earth. 


RESEARCH  ACCOMPLISHED 


Our  NACT  shear  wave  model  (Figure  1)  has  been  developed  from  a  global  data  set  of  surface  wave  waveforms 
crossing  the  region  of  interest  (longitude  30  to  150  degrees  and  latitude  -10  to  60  degrees).  This  region  is  highly 
heterogeneous,  presenting  a  challenge  for  calibration,  but  it  is  well  surrounded  by  earthquake  sources,  and  a 
significant  number  of  high  quality  broadband  digital  stations  exist.  We  started  from  the  existing  waveform  database 
that  was  collected  at  Berkeley  over  the  last  10  years  for  the  construction  of  global  mantle  tomographic  models  (Li 
and  Romanowicz,  1996;  Megnin  and  Romanowicz,  2000;  Gung  et  al.,  2003;  Panning  and  Romanowicz,  2004)  and 
added  data  from  -20  new  events  in  the  period  up  to  2005.  We  now  have  38826  3-component  waveforms  from 
393  events  recorded  at  169  stations.  The  data  has  been  processed  using  an  automated  algorithm,  which  removes 
glitches  and  checks  for  many  common  problems  related  to  timing,  poor  instrument  response,  and  excessively  noisy 
windows.  A  weighting  scheme  has  been  applied  to  insure  even  distribution  of  data  across  the  region.  This  model  is 
parameterized  laterally  in  spherical  spline  level  6,  which  correspond  to  lateral  resolution  of  -200  km.  The 
corresponding  radially  anisotropic  model  is  parameterized  in  the  spline  level  5,  which  corresponds  to  a  lateral 
resolution  of  -400  km. 

The  next  step  in  our  study  has  been  to  perform  an  iteration  of  our  3-D  model  using  the  “non-linear”  Bom 
approximation  (N-Bom),  which  is  modified  from  the  more  standard  3-D  “linear”  Bom  approximation  by  including  a 
“Path  Average”  term.  This  term  allows  the  accurate  inclusion  of  accumulated  phase  shifts  that  arise  in  the  case  when 
the  wavepath  crosses  a  spatially  extended  region  with  a  smooth  anomaly  of  constant  sign.  The  linear  Bom  terms 
account  for  single  scattering  effects  outside  of  the  great  circle  path  and  are  modeled  according  to  the  expressions  of 
Capdeville  (2005).  Accounting  for  scattering  outside  of  the  great-circle  path  is  the  one  difference  with  our  initial 
NACT  approach.  We  demonstrate  the  importance  of  the  non-linear  effect  in  Figure  2  by  comparing  3-D  Bom  and 
3-D  N-Bom  synthetics  with  those  of  SEM.  This  particular  case  illustrates  the  importance  of  including  the  non-linear 
“path  average”  phase  shift  for  paths  that  accumulate  large  phase  delays. 

Because  the  calculation  of  the  3-D  Born  sensitivity  kernels  is  very  expensive  computationally,  we  have  to  select  a 
subregion  (longitude  75  to  150  degrees  and  latitude  0  to  45  degrees),  where  the  lateral  heterogeneities  are  more 
significant  than  the  surrounding  regions.  In  order  to  further  reduce  the  intensity  of  the  computation,  we  require  that 
both  events  and  stations  must  be  in  the  large  region  (longitude  30  to  150°  E  and  latitude  -10  to  60°  N)  and  only  the 
ray  paths  along  the  minor  arcs  are  selected.  We  calculated  3-D  Born  sensitivity  kernels  for  162  events  using  the 
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computing  facilities  (Jacquard)  of  the  National  Energy  Research  Scientific  Computing  Center 
(NERSC — www.nersc.gov).  When  we  generate  synthetics  for  the  N-Bom  inversion,  we  use  PAVA 
(Woodhouse  and  Dziewonski,  1984)  in  the  subregion  to  obtain  the  non-linear  phase  term  and  NACT  outside  the 
subregion.  The  adopted  damping  scheme  for  isotropic  and  radially  anisotropic  models  is  the  same  as  that  used  in  the 
NACT  inversion. 


d  InVC^ 


Figure  1.  Starting  isotropic  shear  velocity  model  derived  using  NACT,  shown  here  only  in  the  subregion  (Ion 
75  to  150  degrees  and  lat  0  to  45  degrees)  with  map  grid  interval  of  15  degrees.  The  lateral 
parameterization  is  in  level  6  spherical  splines,  which  give  a  resolution  of  ~200  km.  Values  are 
shown  as  percent  perturbations  from  the  isotropic  velocity  of  the  Preliminary  Reference  Earth 
Model  (PREM). 


Our  N-Bom  shear  velocity  model  is  shown  in  Figure  3.  Both  N-Bom  and  NACT  derived  models  can  fit  waveforms 
very  well,  with  up  to  -83%  variance  reduction  (depending  on  the  choice  of  damping),  and  the  N-Bom  residual 
variance  is  even  lower  than  that  of  NACT  (Figure  4).  While  the  models  agree  in  general,  there  are  some  notable 
differences  between  them  in  detail.  For  example,  beneath  the  Tibetan  plateau,  the  N-Bom  model  shows  a  stronger 
fast  velocity  anomaly  in  the  depth  range  150  km  to  250  km,  which  disappears  at  greater  depth.  This  indicates  that 
there  is  no  delamination  of  lithosphere  beneath  the  plateau,  as  has  been  suggested  by  some  authors. 
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■BORN 

■NBORN 


Figure  2.  A  comparison  of  synthetic  seismograms  using  different  approximations  for  paths  through  a  large 
slow  anomaly.  The  paths  (numbered  1-3)  projected  on  a  slice  of  the  model  at  220  km  depth  are 
shown  at  top.  The  top  row  shows  SEM  synthetics  through  the  1-D  background  model  (SEM1D; 
dashed)  and  through  the  3-D  model  (SEM3D;  solid)  paths  1  (left),  2  (middle),  and  3  (right).  The 
middle  row  shows  the  predicted  3-D  synthetics  plotted  on  top  of  the  SEM3D  trace  (black)  for  the 
path-average  approximation  (PAVA;  red),  linear  Born  approximation  (BORN;  green),  and  non- 
linearly  modified  Born  (NBORN;  blue).  The  bottom  row  shows  the  misfits  for  PAVA,  BORN,  and 
NBORN  to  the  SEM3D  trace.  The  amplitude  is  scaled  to  the  maximum  of  the  SEM3D  trace  in  each 
column.  While  the  Born  approach  is  important  to  represent  the  3-D  effects  in  many  configurations, 
this  figure  illustrates  the  case  of  a  path  sampling  the  center  of  a  large  anomaly.  In  this  case,  the  non¬ 
linear,  path  average  contribution  to  the  seismogram  dominates,  and  the  linear  Born  approximation 
does  a  poor  job  at  predicting  the  waveforms. 
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Figure  3.  N-Born  isotropic  shear  velocity  model  derived  using  the  “non-linear”  3-D  Born  approximation  in 
the  same  subregion  as  in  Figure  1. 


Figure  4.  Residual  variance  as  a  function  of  the  damping  factor.  In  both  N-Born  and  NACT  models,  the 
chosen  damping  factor  for  the  models  shown  respectively  in  Figure  3  and  Figure  1  is  5.0E-5. 

In  collaboration  with  researchers  at  Institut  de  Physique  du  Globe  de  Paris  (IPGP),  we  have  completed  a  regional 
version  of  the  Spectral  Element  code  (RegSEM.l).  The  regional  SEM  incorporates  PML  boundary  conditions  on  the 
lateral  borders  of  the  region,  which  effectively  eliminate  spurious  reflections  from  the  borders,  and  a  non- 
conforming  mesh.  Compared  with  previous  versions,  first  of  all,  both  Moho  topography  and  surface  topography 
have  been  included  into  the  process  of  building  the  mesh.  With  reference  to  a  Moho  model,  every  element  in  the 
simulation  domain  is  assigned  an  integer  (1  for  elements  just  above  the  Moho,  -1  for  those  just  below,  and  0  for  the 
others),  and  we  are  able  to  provide  realistic  velocity  contrasts  at  the  Moho  discontinuity  while  introducing  any  3-D 
elastic  mantle  model.  Second,  anisotropy  has  been  included  in  this  new  version  of  regional  SEM  (any  anisotropic 
structure  can  be  considered).  Third,  we  do  not  need  to  rotate  events,  stations,  and  models  to  be  centered  at  the  North 
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Pole  before  the  simulation,  as  required  in  the  previous  versions.  All  the  rotations  can  be  done  automatically,  and  all 
the  input  information  related  to  the  event  focal  mechanism  and  station  coordinates  are  straightforward.  Thus,  it  is 
now  very  convenient  to  add  more  stations  into  the  station  list.  Finally,  this  new  version  of  regional  SEM  has  been 
fully  optimized.  The  same  computation  now  takes  only  one-third  of  the  previous  cluster  hours.  This  will  be  very 
helpful  for  our  future  waveform  inversion  using  the  regional  SEM  method.  For  demonstration,  we  examine  the 
effect  of  the  Moho  topography  (CRUST2.0)  on  the  synthetic  waveforms  computed  in  our  NACT  velocity  model 
(Figure  5).  Inclusion  of  the  Moho  topography  can  significantly  change  the  fundamental  Rayleigh  waveforms  in  not 
only  travel  time  but  also  amplitude.  For  the  ray  paths  with  significant  oceanic  portion  (DAV),  the  shape  of  the 
waveform  is  also  changed  (Figure  6). 


Figure  5.  Map  showing  event  and  station  locations  for  the  SEM  simulation  shown  in  Figure  6.  The 
background  is  our  most  recent  model  derived  using  NACT,  shown  at  a  depth  of  60  km. 


Figure  6.  Comparison  of  vertical  component  SEM  synthetics  in  our  NACT  derived  model  (Figure  1)  with 
(red)  and  without  (black)  consideration  of  the  Moho  topography.  Synthetics  have  been  low-pass 
filtered  with  a  corner  at  60  sec. 
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In  order  to  refine  the  velocity  structure  beneath  the  region  of  our  study,  we  perform  forward  waveform  modeling 
with  the  FKI  method.  Broadband  seismograms  are  downloaded  from  the  Incorporated  Research  Institutions  for 
Seismology  (IRIS)  and  corrected  to  absolute  ground  velocity  (cm/sec).  We  show  2  event  locations  of  events 
2000256  (9/12/2000,  Mw6.1)  and  2003107  (4/17/2003,  Mw6.3)  and  the  IRIS  station  distributions  (Figure  7).  We 
start  with  the  2000256  event,  for  which  the  continental  ray  paths  are  dominant,  to  obtain  the  1-D  velocity  structure 
between  the  source  and  each  receiver.  Broadband  data  are  bandpass  filtered  at  0.005-0.05  Hz.  We  used  the  Harvard 
centroid  moment  tensor  (CMT)  solution  for  the  source  parameters,  and  the  starting  model  is  a  1-D  layered  average 
crustal  velocity  structure  derived  from  CRUST2.0.  Using  the  best  velocity  model  we  can  obtain  (Figure  8),  we 
compute  Green’s  functions  and  perform  the  moment  tensor  analysis  for  two  ranges  of  frequency  (0.01-0.05  Hz  and 
0.005-0.03  Hz).  Then,  we  select  the  event  2003107,  for  which  we  have  ray  paths  similar  to  event  2000256,  to 
perform  the  moment  tensor  analysis  using  Green’s  function  obtained  from  our  1-D  simulation  (Figure  9).  We  find  a 
moment  tensor  solution  in  good  agreement  with  the  CMT  solution,  whereas  the  solution  obtained  using  the  PREM 
reference  model  is  very  poor.  While  this  example  was  chosen  because  we  expect  that  we  can  use  Harvard  CMT 
solutions  for  M>6  events  as  good  references,  this  indicates  that  the  additional  regional  modeling  effort  is  worthwhile 
and  will  lead  to  better  moment  tensor  solutions  for  smaller  events  in  the  area,  when  we  extend  the  modeling  to 
higher  frequencies  (0.02-0.05  Hz). 


90°  120° 


Figure  7.  Event  2000256  (9/12/2000,  Mw6.1),  2003107  (4/17/2003,  Mw6.3),  and  IRIS  station  distributions. 
Source  2000256  (9/12/2000,  Mw6.1)  is  used  to  obtain  the  velocity  structures. 
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Figure  8.  Best  P-wave  and  S-wave  velocity  structures  for  the  paths  between  event  2000256  and  IRIS  stations 
obtained  from  1-D  forward  modeling. 
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2003107 


Figure  9.  Top:  Moment  tensor  solution  for  event  2003107  in  the  0.005-0.03  Hz  range  using  Green’s  function 
computed  from  the  velocity  model  obtained  using  the  2000256  event  (Figure  9).  Bottom:  Global 
CMT  solution  from  Harvard  University. 
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CONCLUSIONS  AND  RECOMMENDATIONS 

A  3-D  non-Linear  Born  implementation  is  shown  to  be  promising  to  accurately  model  large  scale  lateral  structures. 
This  approach  includes  relatively  rapid  calculation  of  partial  derivatives  for  use  in  the  next  inversion  step  where  the 
forward  modeling  will  be  performed  using  the  recently  completed  regional  SEM  numerical  code. 

Further  work  on  the  regional  SEM  synthetics  and  3-D  Born  sensitivity  kernels  should  offer  continued  improvement 
of  the  model.  We  are  in  the  process  of  implementing  an  iteration  of  our  model  using  the  regional  SEM  code. 
Additionally,  other  approaches  and  data  sets,  including  ground-truth  data  sets  from  the  Knowledge  Base,  teleseismic 
receiver  functions,  broadband  waveform  forward  modeling  at  shorter  periods,  and  surface  wave  group  velocities  will 
allow  for  validation  and  improvement  of  the  model  and,  in  particular,  for  improved  moment  tensor  solutions. 
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